Separation of two regimes in a disordered boson Hubbard Model 
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We study the insulator-to-superfluid transition in a two-dimensional disordered boson Hubbard 
model at zero temperature for intermediate strength of disorder at commensurate density. Via Monte 
Carlo calculations of the correlation functions in the integer current representation of the model, 
we obtain the dynamical critical exponent z = 1.5 ± 0.1, supporting the multi-critical behavior 
separating the strong and weak disorder regimes. Investigating the density profile, we suggest that 
the density fluctuations due to the particle-hole excitations drive the transition in the weak disorder 
regime. 
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Localization in strongly correlated systems has at- 
tracted attentions as an important problem in condensed 
matter physics. The interplay of disorder and interac- 
tion in these systems still remains not fully understood. 
Among those models adopted to investigate the role of 
disorder, the boson Hubbard model is one of the most 
convenient models partly because the order parameter is 
well defined. The insulator-to-superfluid transitions of 
this model 1] in low dimensions due to both interaction 
and disorder have been studied intensively in the past 

decade HHSHE0EHII3- ThiS problem also has 
many experimental counterparts such as 4 He in p orous 

media [ll| and cold atoms in optical lattices |12U13j . Also 

the superconducting transitions at zero temperature in 

thin films [lj] are believed to belong to the universality 

classes represented by the boson Hubbard model. The 

phase transitions in small Josephson- j unction arrays [lfij 

and granular superconductors [lfj have been described by 

similar models. 

The boson Hubbard model can be written as 
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where b r (bl) is the boson destruction (creation) operator 
at site r, n r = b\b r is the boson number operator, and 
U and t are the on-site repulsion energy and the hop- 
ping matrix element, respectively. Here we assume that 
disorder is given in the local chemical potential so that 
jl — A < /i r < p, + A, where jl is the average chemical 
potential which controls the total density of bosons. The 
parameter A characterizes the strength of disorder. Dif- 
ferent realizations of disorder in random magnitude 
or random phase of the hopping matrix t have been 
discussed recently. 

In pure case (A = 0), the insulating phase of the model 
has a finite Mott energy gap for excitations, and the den- 
sity fluctuations are suppressed. The Mott gap vanishes 
at the transition. The relevant density fluctuations driv- 
ing the phase transition, as t increases, for the commen- 
surate density (jl — 0) are the particle- hole excitations. 



This nature of the transition results in the dynamical 
critical exponent 2=1. For an incommensurate density 
(p 7^ 0), however, either single particle or hole excita- 
tions are favored since the chemical potential breaks the 
particle-hole symmetry globally. In this case, we have 
2 = 20. 

In the presence of disorder (A ^ 0), the situation is 
subtle. Since non-zero /i r breaks the particle-hole sym- 
metry locally, by adopting a single-particle picture, it has 
been argued that an arbitrarily weak random potential 
will localize a single-particle excitation as soon as the 
Mott gap vanishes. In other words, near the transition 
where the Mott gap is smaller than A, it is always pos- 
sible to find an arbitrarily large region in which the 
chemical potential is uniformly shifted to make the single- 
particle excitation gap vanish. Therefore, there exists an 
insulating phase, called Bose glass (BG) phase, which has 
finite regions with vanishing gap. Phase transition occurs 
when this region grows infinitely. A hallmark of this tran- 
sition is the fact that z = d (d is the dimensionality of 
the system), since the BG phase is compressible 0. 

This scenario of the BG-to-SF transition governed by 
growing zero-gap regions is certainly convincing for large 
jl or large A since particles or holes cause the density 
fluctuations separately. However, near jl — for weak 
disorder, there is a possibility that the particle-hole exci- 
tations still survive to drive the transition possiblely due 
to the existence of the so-called statistical particle-hole 
svmmetrv|20j. i.e., that this symmetry is restored when 
a random potential is averaged on the scales of diverging 
correlation length. Therefore, it is an interesting ques- 
tion to check whether the particle-hole excitations govern 
the transition in the presence of weak disorder to change 
the nature of the transition from that of the BG-to-SF 
transition. 

In the quantum rotor model equivalent to the boson 
Hubbard model in the limit of large density, a direct tran- 
sition from the Mott insulator (MI) to SF transition was 
reported[2ll for a weak disorder A = 0.2 at jl = with 
2=1, instead of the BG-to-SF transition. Subsequently 
some numerical evidences were provided, which supports 
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that the direct MI-to-SF transition occurs even for non- 
zero small \x l^l- Further Monte Carlo studies}^ of the 
same model suggested that the nature of the transition is 
indeed divided by a multi-critical line into the strong dis- 
order regime where the BG-to-SF transition occurs and 
the weak disorder regime where the direct MI-to-SF tran- 
sition occurs. As an evidence for the existence of the 
multi-critical line which defines the critical strength of 
the disorder, A c (ft), the new value z — 1.35 ± 0.05 was 
found on the line through the finite-size scaling behav- 
ior of the superfluid stiffness. Also the possibility that 
weak disorder is screened out by the proliferation of the 
particle-hole excitations near the transition around /2 = 
was raised. 

One of the key quantities closely related with the na- 
ture of the transition is the dynamical critical exponent, 
z. In this work, we calculate z directly from the corre- 
lation functions through Monte Carlo simulations. The 
motivation is that the numerical value of z obtained in 
this way is insensitive to some scaling parameters such as 
the aspect ratio and K as far as K « K c , in contrast to 
the finite-size scaling behavior which is in general very 
sensitive to those parameters. Therefore, it is another 
method to obtain the critical exponents even though its 
accuracy would be poorer than that of a finite-size scal- 
ing method. We also investigate the density profile to see 
how the density fluctuations are raised. 

In the limit of large density, the boson Hubbard model 
can be reduced to a quantum rotor model 
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where J = 2not with the average number of bosons per 
site no, and 9 T is the phase angle of the local order pa- 
rameter. This model has the global U(l) symmetry (i.e. 
the Hamiltonian is invariant under the transformation 
9 r — » 9 r + (j) for a constant </>). It implies that the ground 
state of this Hamiltonian has two phases: one with this 
symmetry, which is an insulating phase, and the other, 
the superfluid (SF) phase, in which the symmetry is spon- 
taneously broken. 

This transition has been studied in an equivalent clas- 
sical model 10] whose partition function, represented by 
integer current vectors J, is given by 
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where x = (r, r) is (d + 1) dimensional lattice point, 
and the expression, V • J — 0, denotes the current con- 
servation condition. Here K ~ \/t/U is the parameter 
controlling the quantum fluctuations. 

The correlation function, in terms of 9 in Equation^, 
is defined as 

C(ri,r 2 ,Ti,r 2 ) = [(cxp{i6»(ri, n) - i9(r 2 , T 2 )})] a v, (4) 



where (• ■ •} denotes the ensemble average and [• • -] av de- 
notes the average over different disorder complexions 
with the same strength. In terms of the current vari- 
ables, it can be rewritten as 
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for the equal-time correlation functions, and 
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for the zero-distance correlation functions. (r) are the 
correlation functions for t > or t < 0. When fi = 0, we 
expect C^r (r) = C T (r) , implying the global particle- hole 
symmetry. Here a path is chosen which connects two 
points separated in distance r in a spatial direction or t 
in the temporal direction. Recently a new algorithm was 
developed[3 E3 f° r calculating these correlation func- 
tions very efficiently. By adopting this algorithm, the 
correlation functions are calculated via Monte Carlo sim- 
ulations in (2+1) dimensional lattices near K c w 0.292 
when p, = and A = 0.4. This is the critical strength of 
disorder within the error range. The critical point K c is 
determined from the finite-size scaling of the superfluid 
stiffness. Typically the average is taken over 1,000 - 2,000 
disorder realizations. For each disorder configuration, we 
construct 2,000 worms and 40,000 worms for the smallest 
system (8 x 8 x 17) and the largest system (24 x 24 x 73) 
in our simulations, respectively. 

The results are shown in Figures ^ an d [21 m different 
sizes and for different K's. In order to extract the critical 
exponents, we fit the data on scaling curves. The scaling 
form for the correlation function is given by 

C(r, r) = (r 2 + ar 2 /*)-( d +*- 2 +"V2 /(r / Ci T / e ), (7) 

where / is a scaling function and o is a non-universal 
constant. At the critical point (£ — > oo), the correlation 
functions in finite size systems will have forms 
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where c and d are some constants, £ x and £ T are param- 
eters characterizing the correlation lengths in x and r 
directions respectively, y x = d+z — 2 + rj, and y T = y x /z. 
By using these forms, we obtain y x = 1.44 ± 0.03 and 
y T = 0.92±0.07 at if = 0.292, implying that z= 1.5±0.1. 
£e grows rapidly as K — > 0.292 from below whereas 
£ r ~3> L T for K shown in the figures. The value of 77 
is too uncertain to estimate a reliable number. The error 
ranges of y x and y T are estimated as the range of best 
fittings for different sizes and aspect ratios. 



This value of the dynamical critical exponent is nearly 
consistent with the previously reported value z = 1.35 ± 
0.05. Note that we expect z = 2 for strongly disordered 
cases and z = 1 for weakly disordered cases, if disor- 
der is irrelevant, in two dimensions as discussed above. 
Therefore, this anomalous value of z suggests that A c de- 
termines the multi-critical line separating the weak and 
strong disorder regimes. 

What makes these regimes different is an interesting 
question. Here we propose that the relevance of the 
particle-hole excitations to the transition is different in 
these two regimes. In order to check this conjecture, we 
plot the density profiles. The net particle density at r on 
xy-plane for a given disorder complexion is 



(10) 



We define the number of particles and holes as N p — 
^2 T n r 9(n r ) and Nh = — ^2 T n r 9(—n T ), respectively, 
where 9 is the step function. 

Figure|31is a typical density plot of particles (filled cir- 
cles) and holes (open circles) in the presence of disorder 
for A = 0.4 at K — 0.320. Their density at a certain 
site is denoted by the area of the circles. Similar features 
appear for different A and K, In the Mott insulating 
regime neither particles nor holes appear. The fluctua- 
tions of particles and holes grow as K increases. However, 
they show different behavior in weak and strong disorder 
regimes, which is summarized in Figure 0] 

Figure0]shows the number of particles and holes when 
p, = for A = 0.3, 0.4, and 0.5 in a 10 x 10 x 25 lattice. 
The same complexion of disorder is used with different 
A. Note that the critical points in these systems are 
K c w 0.313,0.292, and 0.248 for A = 0.3,0.4, and 0.5, 
respectively, and the critical strength of the disorder is 
A c w 0.4. We see that, for weak disorder (A = 0.3), 
the number of particles and holes matches each other. It 
strongly suggest that in this regime the low-energy ex- 
citations are the particle-hole pairs. In addition, their 
fluctuations increase sharply around the transition, im- 
plying that a finite compressibility rises at the transi- 
tion as the energy for density fluctuations vanishes. All 
these features supports that a direct MI-to-SF transi- 
tion actually happens in this regime. For strong disorder 
(A = 0.5), on the other hand, the numbers of particles 
and holes change separately. It means that the single 
particle or hole states are occupied separately. The tran- 
sition occurs as delocalized single particle states are oc- 
cupied. Also the fluctuation increases linearly near the 
transition. This implies that there is no abrupt change 
of the compressibility at the transition, consistent with 
the compressible BG insulating picture. The occurrence 
of the BG-to-SF transition in this regime has been con- 
firmed [Icj . (Note that whc n A = 0.5, K c is independent 
of p.) Therefore, the different behavior of particle and 
hole fluctuations as a function of A as shown in Figure 01 
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FIG. 1: The correlation functions on L x L x L T = 24 x 24 x 73 
lattice near K — 0.292. By fitting the data, we find that 
y x = 1.43 and y T = 0.97 at K = 0.292. 



supports the existence of the critical strength of disorder 
dividing the weak and strong disorder regimes. 

In summary, we calculate the correlation functions of 
the disordered boson Hubbard model in the current rep- 
resentation in (2+1) dimensions via Monte Carlo simu- 
lations. For intermediate strength of disorder (A = 0.4) 
with commensurate filling (fi = 0), we obtain the dy- 
namical critical exponent z = 1.5 ± 0.1, consistent with 
the finite-size scaling behavior of the superfluid stiffness. 
This supports the existence of the multi-critical line sep- 
arating the strong and weak disorder regimes. The num- 
ber of particles and holes as a function of K strongly 
supports the scenario that the density fluctuations due to 
the particle-hole excitations, rather than a single-particle 
excitations, indeed drive the phase transition in weak dis- 
order regime. On the contrary, in the strong disorder 
regime, single particle (hole) states are relevant to the 
transition. 

This work was supported by Korea Research Founda- 
tion Grant (KRF-2002-041-C00110). 
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FIG. 3: A density plot for commensurate density at K — 
0.320 and A = 0.4. Filled (open) circles denote particles 
(holes). Their density at a certain site is proportional to the 
area of the circle centered on that site. There are narticles and 
holes, and 1 reases. 
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FIG. 4: Net number of particles (N p : filled triangles) and 
holes (Nh- open circles) in a 10 x 10 x 25 lattice as a function 
of K. The same complexion is used with different A. Lines 
are guide for eyes. 
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